Flux distribution of metabolic networks close to optimal biomass production 
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We study a statistical model describing the steady state distribution of the fluxes in a metabolic 
network. The resulting model on continuous variables can be solved by the cavity method. In 
particular analytical tractability is possible solving the cavity equation over an ensemble of networks 
with the same degree distribution of the real metabolic network. The flux distribution that optimizes 
production of biomass has a fat tail with a power-law exponent independent on the structural 
properties of the underling network. These results are in complete agreement with the Flux-Balance- 
Analysis outcome of the same system and in qualitative agreement with the experimental results. 



Recently large attention has been addressed by the 
physics community to critical phenomena [l[ in complex 
networks 0, 0, EL H @, . 

The complex topologies, usu- 
ally characterized by non-Poisson degree distributions, 
have large effect on the critical point and critical expo- 
nents of the dynamical models defined on them. The 
Ising model d, H, the epidemic spreading [ll[ and 
the synchronization dynamics [l2| are examples of dy- 
namics models, where the complex structure has strong 
implications. Furthermore in the last decade we have 
assisted to a big breakthrough in system biology, the in- 
terdisciplinary field that studies the biological problems 
going beyond the single biomolecule framework, with the 
description of the intertwined reactions between the con- 
stituents of the cell in terms of networks. This has gener- 
ated a new theoretical framework in which new biological 
statistical findings have been formulated [HI, [3 H5j . In 
system biology there was also the fast development of 
"in silico" biology in which new experiments are simu- 
lated and the predictions are made to stimulate further 
experimental confirmations of the phenomena. A key 
example of a biological system in which the network pic- 
ture is crucial and the "in silico" biology has made rel- 
evant advances, is the prediction of the growth rate of 
single cells of different organisms and the study of the 
metabolic networks. Two key advances in this field have 
been the full characterization of the chemical reactions 
[ItJ for a series of model organisms, as different strain 
of Escherichia coli and Saccaromyces cerevisiae (see for 
example the BIGG database [IB]), and the application 
of the techniques of linear programming for the study of 
the flux of the reaction, an extension which goes under 
the name of Flux-Balance- Analysis [17[ • 

The set of stoichiometric interactions in the cell can be 
represented as a network whose nodes are of two types, 
the metabolic substrates of metabolites and the nodes 
representing the reactions. This bipartite network goes 
under the name of factor graph. In a factor graph, to 
each reaction i is assigned a flux variable s< and to each 
metabolite /i is assigned a steady state condition for the 
production/consumption of the metabolites. The struc- 
ture of the metabolic network has a projection on the 
metabolites which has a power-law degree distribution 
[lH and a hierarchical structure 0,ElJ- I n the metabolic 
networks, to each reaction corresponds an enzyme which 



regulates the rate of each reaction and modulates the flux 
of the reactions. Consequently the maximal flux rate is 
fixed by the maximal enzyme concentration inside the 
cell. Solving the non-linear mass-law equations is a hard 
problem in networks of thousand of nodes. To overcome 
this problem in Flux-Balance- Analysis for each reaction 
a new variable, its flux, is introduced. Each flux includes 
all the dynamical effects associated to each reaction of 
the organism. The Flux-Balance-Analysis [TtI H3] con- 
siders the steady state of the dynamics which optimizes 
the production of the biomass by linear programming. 

The underline assumption of Flux-Balance- Analysis, 
that the cell organisms optimize the biomass is very well 
confirmed by experimental results (20| conducted in rich 
medium but is not fully supported by experiments mea- 
suring the growth rate of single knockout strains. For 
this reason other algorithms have been designed relaxing 
this condition. [2III22I]. 

In this paper we will study the flux distribution in 
metabolic networks that has a heavy tail as found in 
experiments [23| and in Flux-Balance- Analysis [24| pre- 
dictions in Echerichia coli. In particular we add to the 
description of the metabolic networks some theoretical 
statistical mechanics insights using the cavity method 
0, HB]. Different theoretical models have been already 
proposed [H, H3j HI] for the flux distribution but neither 
of them has been able to theoretically predict the out- 
come of the experiments or of the Flux-Balance- Analysis 
calculations. Here we will relate the power-law exponent 
of the flux distribution with the steady distribution as 
an indication of criticality. In fact it can be predicted by 
assuming that the network is optimizing the biomass pro- 
duction. Thus we will characterize this optimized state 
doing an asymptotic expansion of the cavity equation 
close to the critical point and we will measure the crit- 
ical exponents corresponding to this critical transition. 
The method which we formulate here is a new method to 
solve the cavity equation on continuous variables defined 
on a compact interval of the real axis and it can be ex- 
tended to other critical phenomena on complex networks 
and continuous variables defined in a limited interval. We 
find that the distribution of the fluxes present in the op- 
timized state develops a power-law tail with exponents 
that are independent on the structure of the underlying 
network, a new scenario in which the complex topology 
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FIG. 1: The cavity tree C^. If the metabolite \i is removed 
from the metabolic factor graph, to each reaction i in which it 
takes part we can assign a flux st with the cavity distribution 




FIG. 2: The degree distribution of p(k) and p(q) for _Bs- 
cherichia coli, data taken form the BIGG database [l6l].The 
line indicates the power-law p(k) = k ' with 7 = 3.0. 



does not affect the critical state of the cell. The power- 
law exponent that we find is in full agreement with the 
results of Flux-Balance- Analysis [24| and only partially in 
agreement with the experimental results [231 ] 

The model- The metabolic network has a bow tie struc- 
ture (l9| : the metabolites are divided into input metabo- 
lites which are provided by the environment, the out- 
put metabolites which provide the biomass and inter- 
mediate metabolites. The stoichiometric matrix is given 
by ((£f )) where [i — 1, . . . , M indicates the number of 
metabolites and i — 1,...,N the number of reactions 
and the sign of & indicates if the metabolite \i is an in- 
put or output metabolite of the reaction i. In the Flux- 
Balance- Analysis method we assume that each interme- 
diate metabolite has a concentration c M at steady state, 
i.e. 



dt / " 



= 



(1) 



where Sj is the flux of the metabolic reaction i. For the 
metabolites present in the environment and the metabo- 
lites giving rise to the biomass production we can fix the 
incoming flux given by 



dc^ = y 



in /out 

9 u 



(2) 



We have already mentioned that the fluxes have some 
biological limitations. To describe these limitations we 
assume that the fluxes G [0, L}. 

The volume of solutions V of this problem is given by 



V 



/ ... / Y[d Sl Y[S(J2^^ s i- 9^)- 
•f «/ „• 1 , . 



(3) 



Belief Propagation- Let us assume that the factor 
graph of the metabolic network has a local tree like struc- 
ture as shown in figure [TJ In this assumption Belief 



Propagation (BP) equations are able to fix the proba- 
bility distribution of the metabolic fluxes with the mea- 
sure defined in ([3]). BP equations are defined on cavity 
graphs. The cavity graphs C M is the full factor graph of 
the metabolic network in the absence of metabolite /1. In 
Cfj, the flux Si of a reaction i in which p is reacting has 
a cavity distribution pi^^(si). Expressing Pi^^{si) in 
terms of the cavity distribution pj^ v {si) where i,j,p,,v 
are defined in figure [TJ we get the BP equations: 



Pi->p(Si) 
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Solving the BP equations for the cavity distributions 
the marginal probability of a flux s, is given by 



Pi(Si) 



a 
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The distribution of the fluxes producing/consuming the 
metabolite p , i.e. £„ = {si}i^N(fi) is given by 



c, 



— s (52tj,» s J - 9») n pj^( s s) ( 5 ) 
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The entropy of the metabolic network can be expressed 
as 



y^(fcj - 1) J dsiPi(si)logpi(si). (6) 
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The optimized state of a random metabolic network- 
We assume that the metabolic network is a random graph 
with M metabolites with degree distribution p(k) and N 
reaction nodes with degree distribution p(q) . In this net- 
work the total number of links is given by N(q) = M(k). 
In Figure [2] we show an example of p{k) and p(q) distri- 
butions for Escherichia coli. The p(k) degree distribution 
for this organism is power-law p{k) ~ k~ J with an ex- 
ponent 7 ~ 3 while the p(q) distribution is much more 
picked. In different organisms the distribution of p(q) 
and p(k) do change but the general scenario of power- 
law p{k) distribution and finite scale p{q) distribution 
remains unchanged. 

To solve the BP equations (j4|) in this random ensemble 
of graphs with given p{k) and p{q) degrees distributions 
we introduce the quantity 



P(s) 



{q)N 



.(*) 



P(g) 



(7) 



dom sign of the reactions and the outgoing and ingo- 
ing flux distribution i.e. over the distribution for the 

Ul P({0) = IL FU* § [Sfou - 1) + S^i,u + 1)] and 
(...) indicates the average over the probability distribu- 
tion p(g) of the g's. In particular we take 

p(g) =PiS(g + gi) +P2S(g-g 2 ) + (1 -pi -P2)S(g) (8) 

where p\ and pi indicate respectively the fraction of in- 
put metabolites and the fraction of metabolites in the 
biomass definition and where <?i is fixed by the environ- 
ment conditions and 172 is the rate of the biomass. 

The Fourier transform of P(s) is the function x( w )-> 



X(w) = / dse iws P(s). 



where the overline m indicates the average over ran- consisted equation 
I 



(9) 



where w = =f~n. The function x( w ) satisfies the self 
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*(X>„(-1) 



p(g) 



w). 



where (•} indicates the average over the distribution of 
the incoming/outgoing fluxes g and where and we have 
assumed I » 1 allowing for a continuous variables w„, 
w. 



C({k v }) 




1 K 



(10) 



The equation for x( w ) wn l have solutions until the rate of 
biomass production g% < G c with G c been correspond- 
ing to the maximal allowed biomass production of the 
metabolic network. 

Close to the critical point gi ~ G c we suppose that the 



J 



distribution P(s) will scales as 

P(s) = s- T <f>(s\g 2 - G c \°) 



(11) 



where $ is a scaling function. In the limit of large L we 
can assume that w (together with the w v ) is a continu- 
ous variable and we can do an asymptotic expansion is 
reflected the scaling behavior of x( w )i i- e - 



x(w) = l-\w\ T - 1 h(w/\g 2 -G c n. 



(12) 



Solving then the self consistent equation for x( w ) Eq. 
(|10[) for analytic distribution p(q) and the distribution of 
the metabolites connectivity decaying like a power-law 
p(k) ~ fc~ 7 . Close to the phase transition we have 



C{{K}) 



9-1 

v=l 



l-{k v - l)|w 1/ |' r - 1 Re h + hku - l)(fc„ - 2)|^| 2(r - 1 )(Re hf 

V 

1 + E IMg^gwKK - l){k v , - 1) + A 2 {g u ,g v .){k v - l)(k u ~ 2){k v , - 1)] 
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- A 3 (g u ,g v ,,g„„)(k v -l)(k v ,-l)(k v „-l) (13) 

v,nu' ,v" 

where Ax, A2 and A3 are linear functions of g v: g v < and g v n , If we develop (jlOp around the point w = 0, we get 



k„ {n„} 



*(5>„(-ir 



w 



(14) 



Since the sum over the connectivity k v are convergent 
in a sparse network, all metabolic networks have for the 
flux distribution the mean field exponents r = 3/2 and 
a = 2 as long as the hypothesis of Flux-Balance- Analysis 
are satisfied. 

The entropy goes like 

X=\g 2 -G c \ a (15) 

with a = er(r — 1) = 1 

Therefore in the physical range for each degree distri- 
bution p(q) or p(k u ) [15| the predicted power-law critical 
exponent for the flux distribution is r = 3/2 in good 
agreement with the Flux-Balance calculations [24| 



Conclusion- We have presented a statistical mechanics 
approach to study the steady state distribution of the 
fluxes in the metabolic network assuming optimization 
of the biomass. The analytic treatment finds a distribu- 
tion of the fluxes which is a power-law with an mean- 
field exponent r = 3/2 independent on the structure of 
the metabolic network. The method can be generalized 
to other critical phenomena defined on continuous vari- 
ables on a finite interval, and work in this direction is in 
progress. 
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